﻿#pragma once

#include "geom_curve.h"
#include "kernel/spline.h"

namespace CAGD
{

/**
 * B 样条曲线.
 */
template <int Dim> class GeomBSplineCurveDim : public GeomCurveDim<Dim>
{
  public:
    /**
     * @brief 次数 + 控制点构造.
     * @param deg
     * @param cpts
     * @param periodic
     * @see
     */
    GeomBSplineCurveDim(int deg, const std::vector<BPntN<Dim>> &cpts, bool periodic = false)
        : m_cpts(cpts), m_periodic(periodic), m_weights(cpts.size(), 1)
    {
        int p = deg;
        int n = cpts.size() - 1;

        // 生成均匀节点向量
        auto &tau = m_knots;
        for (int i = 0; i < p + 1; i++)
            tau.push_back(0);
        for (int i = p + 1; i < n + 1; i++)
            tau.push_back(1.0 * (i - p) / (n + 1 - p));
        for (int i = n + 1; i < n + p + 2; i++)
            tau.push_back(1);
    }

    /**
     * @brief 节点 + 控制点构造.
     * @param knots
     * @param cpts
     * @param periodic
     * @see
     */
    GeomBSplineCurveDim(const std::vector<double> &knots, const std::vector<BPntN<Dim>> &cpts, bool periodic = false)
        : m_knots(knots), m_cpts(cpts), m_weights(cpts.size(), 1)
    {
        // 根据节点和控制点判断周期性
        m_periodic = false;
    }

    /**
     * @brief 获得节点.
     * @return
     * @see
     */
    const std::vector<double> &Knots() const
    {
        return m_knots;
    }

    /**
     * @brief 设置节点.
     * @param knots
     * @see
     */
    void SetKnots(const std::vector<double> &knots)
    {
        m_knots = knots;
    }

    /**
     * @brief 获得节点数量.
     * @return
     * @see
     */
    int NbKnots() const
    {
        return m_knots.size();
    }

    /**
     * @brief 获得控制点.
     * @return
     * @see
     */
    const std::vector<BPntN<Dim>> &Poles() const
    {
        return m_cpts;
    }

    /**
     * @brief 设置控制点.
     * @param cpts
     * @see
     */
    void SetPoles(const std::vector<BPntN<Dim>> &cpts)
    {
        m_cpts = cpts;
    }

    /**
     * @brief 获得控制点数量.
     * @return
     * @see
     */
    int NbPoles() const
    {
        return m_cpts.size();
    }

    /**
     * @brief 次数.
     * @return
     * @see
     */
    int Degree() const
    {
        return m_knots.size() - m_cpts.size() - 1;
    }

    /**
     * @brief 获得周期性.
     * @return
     * @see
     */
    bool Periodic() const
    {
        return m_periodic;
    }

    /**
     * @brief 起始参数.
     * @return
     * @see
     */
    double FirstParameter() const override
    {
        return m_knots[Degree()];
    }

    /**
     * @brief 终止参数.
     * @return
     * @see
     */
    double LastParameter() const override
    {
        return m_knots[m_knots.size() - Degree() - 1];
    }

    /**
     * @brief 取值函数.
     * @param t
     * @param P
     * @see
     */
    void D0(double t, BPntN<Dim> &P) const override
    {
        int p = Degree();
        int l = Spline::FindSpan(t, p, m_knots);

        std::vector<BPntN<Dim>> ps;
        for (int j = l - p; j <= l; j++)
            ps.push_back(m_cpts[j]);

        for (int r = 1; r <= p; r++)
        {
            for (int j = l; j >= l - p + r; j--)
            {
                double lambda = (t - m_knots[j]) / (m_knots[j + p + 1 - r] - m_knots[j]);
                ps[j + p - l] = combine(ps[j + p - l], ps[j + p - l - 1], lambda, 1 - lambda);
            }
        }
        P = ps.back();
    }

    /**
     * @brief 切线函数.
     * @param t
     * @param P
     * @param V1
     * @see
     */
    void D1(double t, BPntN<Dim> &P, BVecN<Dim> &V1) const override
    {
        auto D = DerivPoints(1, t);
        P = D[0];
        V1 = D[1];
    }

    /**
     * @brief 2 阶导函数.
     * @param t
     * @param P
     * @param V1
     * @param V2
     * @see
     */
    void D2(double t, BPntN<Dim> &P, BVecN<Dim> &V1, BVecN<Dim> &V2) const override
    {
        auto D = DerivPoints(2, t);
        P = D[0];
        V1 = D[1];
        V2 = D[2];
    }

    /**
     * @brief n 阶导函数.
     * @param N
     * @param t
     * @param P
     * @param V
     * @see
     */
    void DN(int N, double t, BPntN<Dim> &P, std::vector<BVecN<Dim>> &V) const override
    {
        auto D = DerivPoints(N, t);
        P = D[0];

        V.clear();
        for (int i = 1; i < D.size(); i++)
            V.push_back(D[i]);
    }

    /**
     * @brief 获得引用控制点.
     * @return
     * @see
     */
    std::vector<BPntN<Dim> *> MakeReference() override
    {
        std::vector<BPntN<Dim> *> reference;
        for (int i = 0; i < m_cpts.size(); i++)
            reference.push_back(&m_cpts[i]);
        return reference;
    }

    /**
     * @brief 节点插入算法：将 u 插入 r 次.
     * @param u
     * @param r
     * @see
     */
    void InsertKnot(double u, int r)
    {
        int p = Degree();
        int n = NbPoles() - 1;
        int m = n + p + 1;
        int k = Spline::FindSpan(u, p, m_knots);    // 找到插入位置
        int s = Spline::FindMult(k, u, p, m_knots); // 获得插入节点之前的重数
        std::vector<BPntN<Dim>> &Q = m_cpts;
        std::vector<double> &U = m_knots;

        // 保存用于迭代的控制点 k-p 到 k-s
        std::vector<BPntN<Dim>> R(p + 1);
        for (int i = 0; i < p - s + 1; i++)
            R[i] = Q[k - p + i];

        // 前面 k + 1 个控制点不变，
        // 后面不变的控制点，向后移动 r 位（倒过来移动）
        Q.resize(n + r + 1);
        for (int i = n; i >= k - s; i--)
            Q[i + r] = Q[i];

        // 循环插入
        for (int j = 1; j < r + 1; j++)
        {
            int L = 0; // 用于 R 的索引

            // k - p + r <= i <= k - s
            for (int i = k - p + j; i < k - s + 1; i++, L++)
            {
                double alpha = (u - U[i]) / (U[i + p - j + 1] - U[i]);
                R[L] = R[L + 1] * alpha + R[L] * (1 - alpha);
            }

            // 顶部和底部两个控制点不变
            Q[k - p + j] = R[0];
            Q[k - s + r - j] = R[p - s - j];
        }

        // 存储中间的迭代结果
        int L = k - p + r;
        for (int i = L + 1; i < k - s; i++)
            Q[i] = R[i - L];

        // 最后修改节点
        U.resize(m + r + 1);
        for (int i = m + r; i > k + r; i--)
            U[i] = U[i - r];
        for (int i = 1; i < r + 1; i++)
            U[k + i] = u;
    }

    /**
     * @brief 节点优化算法：将递增序列 X 中的节点插入.
     * @param X
     * @see
     */
    void RefineKnot(const std::vector<double> &X)
    {
        int n = NbPoles() - 1;
        int p = Degree();
        int m = n + p + 1;
        int r = X.size() - 1;
        std::vector<BPntN<Dim>> &Q = m_cpts;
        std::vector<double> &U = m_knots;

        // 获得要插入的区间范围
        double a = Spline::FindSpan(X[0], p, U);
        double b = Spline::FindSpan(X[r], p, U) + 1;

        // 前面 a + 1 个控制点不变
        // 后面的控制点平移 r+1
        Q.resize(n + r + 2);
        for (int j = n; j >= b - 1; j--)
            Q[j + r + 1] = Q[j];

        // 前面 a + 1 个节点不变，后面节点向后平移 r+1
        U.resize(m + r + 2);
        for (int j = m; j >= b + p; j--)
            U[j + r + 1] = U[j];

        // 从最后一个要变化的节点开始向前循环迭代
        int i = b + p - 1;
        int k = b + p + r;
        for (int j = r; j >= 0; j--)
        {
            // 插入一个节点，更新控制点
            while (X[j] <= U[i] && i > a)
            {
                Q[k - p - 1] = Q[i - p - 1];
                U[k] = U[i];

                k--;
                i--;
            }
            Q[k - p - 1] = Q[k - p];
            for (int l = 1; l < p + 1; l++)
            {
                int ind = k - p + l;
                double alpha = U[k + l] - X[j];
                if (fabs(alpha) == 0)
                    Q[ind - 1] = Q[ind];
                else
                {
                    alpha /= U[k + l] - U[i - p + l];
                    Q[ind - 1] = Q[ind - 1] * alpha + Q[ind] * (1 - alpha);
                }
            }
            U[k] = X[j];
            k--;
        }
    }

    /**
     * @brief 计算不大于 d 阶导数的控制点，r1,r2 为索引范围.
     * @return
     * @param d
     * @param r1
     * @param r2
     * @see
     */
    std::vector<std::vector<BPntN<Dim>>> DerivCpts(int d, int r1, int r2) const
    {
        assert(d <= Degree());
        assert(r1 >= 0);
        assert(r2 < int(m_cpts.size()));

        // 0 阶导控制点（初始控制点）
        int r = r2 - r1;
        std::vector<std::vector<BPntN<Dim>>> PK(d + 1, std::vector<BPntN<Dim>>(r + 1));
        for (int i = 0; i < r + 1; i++)
            PK[0][i] = m_cpts[r1 + i];

        int p = Degree();
        for (int k = 1; k < d + 1; k++)
        {
            int tmp = p - k + 1;
            // P_i^k = (P_{i+1}^{k-1}-P_i^{k-1}) / (u_{i+p+1}-u_{i+k}) * (p-k+1)
            for (int i = 0; i < r - k + 1; i++)
                PK[k][i] = (PK[k - 1][i + 1] - PK[k - 1][i]) / (m_knots[r1 + i + p + 1] - m_knots[r1 + i + k]) * tmp;
        }
        return PK;
    }

    /**
     * @brief 计算不大于 d 阶导数曲线在 u 处的点.
     * @return
     * @param d
     * @param u
     * @see
     */
    std::vector<BPntN<Dim>> DerivPoints(int d, double u) const
    {
        int p = Degree();
        std::vector<BPntN<Dim>> CK(d + 1);

        // 导数不超过次数，超过的部分归零
        int du = std::min(p, d);
        int span = Spline::FindSpan(u, p, m_knots);

        // 获得从 0 次开始的所有非零基函数
        std::vector<std::vector<double>> N = Spline::DerivAllBasis(span, u, p, m_knots);
        std::vector<std::vector<BPntN<Dim>>> PK = DerivCpts(du, span - p, span);

        // N[j][r] = N_{r,j}
        for (int k = 0; k < du + 1; k++)
            for (int j = 0; j < p - k + 1; j++)
                CK[k] = CK[k] + PK[k][j] * N[p - k][j];
        return CK;
    }

  protected:
    bool m_periodic;
    std::vector<double> m_knots;
    std::vector<double> m_weights;
    std::vector<BPntN<Dim>> m_cpts;
};

using Geom2dBSplineCurve = GeomBSplineCurveDim<2>;
using GeomBSplineCurve = GeomBSplineCurveDim<3>;

} // namespace CAGD
